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Abstract 

The rotor-router model is a deterministic analogue of random walk. 
It can be used to define a deterministic growth model analogous to 
internal DLA. We prove that the asymptotic shape of this model is a 
Euclidean ball, in a sense which is stronger than our earlier work [T^] . 
For the shape consisting of n = u>df d sites, where ojd is the volume of 
the unit ball in M. d , we show that the inradius of the set of occupied 
sites is at least r — O(logr), while the outradius is at most r + 0(r a ) 
for any a > 1 — l/d. For a related model, the divisible sandpile, we 
show that the domain of occupied sites is a Euclidean ball with error in 
the radius a constant independent of the total mass. For the classical 
abelian sandpile model in two dimensions, with n = nr 2 particles, we 
show that the inradius is at least r/v3, and the outradius is at most 
(r + o(r))/ v / 2. This improves on bounds of Le Borgne and Rossin. 
Similar bounds apply in higher dimensions, improving on bounds of 
Fey and Redig. 

1 Introduction 

Rotor-router walk is a deterministic analogue of random walk, first studied 
by Priezzhev et al. [18] under the name "Eulerian walkers." At each site 
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Figure 1 : Rotor-router aggregate of one million particles in 1? . Each site is 
colored according to the direction of its rotor. 

in the integer lattice Z 2 is a rotor pointing north, south, east or west. A 
particle starts at the origin; during each time step, the rotor at the particle's 
current location is rotated clockwise by 90 degrees, and the particle takes a 
step in the direction of the newly rotated rotor. In rotor-router aggregation, 
introduced by Jim Propp, we start with n particles at the origin; each par- 
ticle in turn performs rotor-router walk until it reaches a site not occupied 
by any other particles. Let A n denote the resulting region of n occupied 
sites. For example, if all rotors initially point north, the sequence will begin 
Ai = {(0,0)}, A 2 = {(0,0),(1,0)}, A 3 = {(0,0), (1,0), (0,-1)}. The region 
^4i,ooo,ooo is pictured in Figure 1. In higher dimensions, the model can be 
defined analogously by repeatedly cycling the rotors through an ordering of 
the 2d cardinal directions in Z rf . 

Jim Propp observed from simulations in two dimensions that the regions 
A n are extraordinarily close to circular, and asked why this was so [7J ITU] . 
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Despite the impressive empirical evidence for circularity, the best result 
known until now [Hi says only that if A n is rescaled to have unit volume, 
the volume of the symmetric difference of A n with a ball of unit volume tends 
to zero as a power of n, as n j oo. The main outline of the argument is 
summarized in [15] . Fey and Redig [5 j also show that A n contains a diamond. 
In particular, these results do not rule out the possibility of "holes" in A n 
far from the boundary or of long tendrils extending far beyond the boundary 
of the ball, provided the volume of these features is negligible compared to 
n. 

Our main result is the following, which rules out the possibility of holes 
far from the boundary or of long tendrils in the rotor-router shape. For 
r > let 

B r = {i£Z d : \x\ < r}. 

Theorem 1.1. Let A n be the region formed by rotor-router aggregation in 
7L d starting from n particles at the origin and any initial rotor state. There 
exist constants c, d depending only on d, such that 

B r —c\ogr C A n C -S r (i_|_ c ' r .- 1 / d logr) 

where r = {n/uJd) l ^ d , and uj^ is the volume of the unit ball in M. d . 

We remark that the same result holds when the rotors evolve according 



to stacks of bounded discrepancy; see the remark following Lemma 5.1 

Internal diffusion limited aggregation ("internal DLA") is an analogous 
growth model defined using random walks instead of rotor-router walks. 
Starting with n particles at the origin, each particle in turn performs simple 
random walk until it reaches an unoccupied site. Lawler, Bramson and 
Griffeath [10] showed that for internal DLA in Z d , the occupied region A n , 
rescaled by a factor of v}' d , converges with probability one to a Euclidean 
ball in M. d as n —> oo. Lawler |11] estimated the rate of convergence. By way 
of comparison with Theorem |1.1| if I n is the internal DLA region formed 
from n particles started at the origin, the best known bounds [11] are (up 
to logarithmic factors) 

B r _ r i/3 C I n C B r+r i/3 

for all sufficiently large n, with probability one. 

We also study another model which is slightly more difficult to define, 
but much easier to analyze. In the divisible sandpile, each site x £ Z d 
starts with a quantity of "mass" uq(x) G M>q. A site topples by keeping up 
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Figure 2: Classical abelian sandpile aggregate of one million particles in Z 2 . 
Colors represent the number of grains at each site. 

to mass 1 for itself, and distributing the excess (if any) equally among its 
neighbors. Thus if x has mass m > 1, then each of the 2d neighboring sites 
gains mass {m — 1) /2d when we topple x, and x is left with mass 1; if m < 1, 
then no mass moves when we topple x. 

Note that individual topplings do not commute; however, the divisible 
sandpile is "abelian" in the following sense. 

Proposition 1.2. Let x\,X2, . . . £ Z rf tea sequence with the property that 
for any x G Z^ there are infinitely many terms x^ = x. Let 

u k (x) = total mass emitted by x after toppling x\, . . . , x^] 
^k( x ) = amount of mass present at x after toppling xi,...,Xk. 

Then u k f u and vj. — > v < 1. Moreover, the limits u and v are independent 
of the sequence {x^}. 

The abelian property can be generalized as follows: after performing 
some topplings, we can add some additional mass and then continue top- 



pling. The resulting limits u and v will be the same as in the case when all 
mass was initially present. For a further generalization, see |16j . 

The limiting function u in Proposition \1.2\ is the odometer function for 
the divisible sandpile. This function plays a central role in our analysis. The 
limit v represents the final mass distribution. Sites x G Z d with v(x) = 1 



are called fully occupied. Proposition 1.2 is proved in section 3, along with 
the following. 

Theorem 1.3. For m > let D m C 1> d be the domain of fully occupied 
sites for the divisible sandpile formed from a pile of mass m at the origin. 
There exist constants c, c' depending only on d, such that 

Br-c C D m C B r+C i , 

where r = {m/u),i) l / d and oj^ is the volume of the unit ball in W 1 . 

The divisible sandpile is similar to the "oil game" studied by Van den 
Heuvel [22]. In the terminology of [5], it also corresponds to the h — * — oo 
limit of the classical abelian sandpile (defined below), that is, the abelian 
sandpile started from the initial condition in which every site has a very 
deep "hole." 

In the classical abelian sandpile model [1] , each site in Z rf has an integer 
number of grains of sand; if a site has at least 2d grains, it topples, sending 
one grain to each neighbor. If n grains of sand are started at the origin in 
Z d , write S n for the set of sites that are visited during the toppling process; 
in particular, although a site may be empty in the final state, we include it 
in S n if it was occupied at any time during the evolution to the final state. 

Until now the best known constraints on the shape of S n in two dimen- 
sions were due to Le Borgne and Rossin [12j . who proved that 

{x G I? | xi + x 2 < vV 12 - 1} C S n C {x G I? | xx, x 2 < Vn/2}. 

Fey and Redig [5] proved analogous bounds in higher dimensions, and ex- 
tended these bounds to arbitrary values of the height parameter h. This 
parameter is discussed in section 4. 

The methods used to prove the near-perfect circularity of the divisible 
sandpile shape in Theorem |1.3| can be used to give constraints on the shape 
of the classical abelian sandpile, improving on the bounds of [5 a and j!2j . 

Theorem 1.4. Let S n be the set of sites that are visited by the classical 
abelian sandpile model in Z rf , starting from n particles at the origin. Write 
n = uodr d . Then for any e > we have 

B Cl r-c 2 C S n C B c ^ r+C ^ 
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Figure 3: Known bounds on the shape of the classical abelian sandpile in 
1? . The inner diamond and outer square are due to Le Borgne and Rossin 
|12j ; the inner and outer circles are those in Theorem 1.4 



where 

ci = (2d-l)- 1/d , c[ = (d-e)- l ' d . 
The constant C2 depends only on d, while c' 2 depends only on d and e. 



Note that Theorem 1.4 does not settle the question of the asymptotic 
shape of S n , and indeed it is not clear from simulations whether the asymp- 
totic shape in two dimensions is a disc or perhaps a polygon (Figure [2| . To 
our knowledge even the existence of an asymptotic shape is not known. 

The rest of the paper is organized as follows. In section 2, we derive the 
basic Green's function estimates that are used in the proofs of Theorems |1.1| 



1.3 and |1.4[ In section 3 we prove Proposition 1 1 .2 1 and Theorem 1.3 for the 
divisible sandpile. In section 4 we adapt the methods of the previous section 
to prove Theorem |1.4| for the classical abelian sandpile model. Section 5 is 



devoted to the proof of Theorem |1.1| 
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2 Basic Estimate 



Write (Xfc)fc>o for simple random walk in Z rf , and for d > 3 denote by 

g(x) = E a #{k\X k = x} 

the expected number of visits to x by simple random walk started at the 
origin. This is the discrete harmonic Green's function in Z rf ; it satisfies 
Ag(x) = for x/o, and Ag(o) = —1, where A is the discrete Laplacian 



&g(x) = wiJ2a{y) -g(x). 



2d 

y~x 

The sum is over the 2d lattice neighbors y of x. In dimension d = 2, simple 
random walk is recurrent, so the expectation defining g{x) is infinite. Here 
we define the potential kernel 

g(x) = lim g n (x) - g n (o) (1) 

n^oo 

where 

g n (x) = E Q #{k < n\X k = x}. 

The limit defining g(x) in ([!]) is finite [HIED], and g(x) has Laplacian Ag(x) = 
for x ^ o, and A<?(o) = —1. Note that ([TJ is the negative of the usual 
definition of the potential kernel; we have chosen this sign convention so 
that g has the same Laplacian in dimension two as in higher dimensions. 
Fix a real number m > and consider the function on Z d 

ld{x) = \x\ 2 +mg(x). (2) 

Let r be such that m = u>dr d , and let 

7d0) = ld{x) -7rf(L r Jei) (3) 

where e\ is the first standard basis vector in Z d . The function 7^ plays a 
central role in our analysis. To see where it comes from, recall the divisible 



sandpile odometer function of Proposition 1.2 



u(x) = total mass emitted from x. 

Let D m C Z d be the domain of fully occupied sites for the divisible sandpile 
formed from a pile of mass m at the origin. For x £ D m , since each neighbor 
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y of x emits an equal amount of mass to each of its 2d neighbors, we have 



Au (x) = ^^2 u (y) -u{x) 

= mass received by x — mass emitted by x 
= 1 - m5 ox . 

Moreover, u = on dD m . By construction, the function 7^ obeys the same 
Laplacian condition: A^fj = 1 — m5 ] and as we will see shortly, 7^ ~ on 
dB r . Since we expect the domain D m to be close to the ball B r , we should 
expect that u ~ 7^. In fact, we will first show that u is close to 7^, and then 
use this to conclude that D m is close to B r . 

We will use the following estimates for the Green's function [61 [2T| ; see 
also [9j Theorems 1.5.4 and 1.6.2]. 

o(x) = /-| lo gM + K + °(M~ 2 )> d = 2 

V ^ \a d \x\ 2 - d + 0{\x\- d ), d>3. K ' 

Here ad = rgz^yj^ , where uod is the volume of the unit ball in M d , and k is 
a constant whose value we will not need to know. For x E Z rf we use \x\ to 
denote the Euclidean norm of x. Here and throughout the paper, constants 
in error terms denoted O(-) depend only on d. 

We will need an estimate for 7^ near the boundary of the ball B r . We 
first consider dimension d = 2. From Q we have 

72(a0 = <P{x) — Km + 0(m\x\~ 2 ), (5) 

where 

<p{x) = \x\ log \x\. 

In the Taylor expansion of 4> about \x\ = r 

4>(x) = <f>(r) - <t>'{r){r - \x\) + ^<f>"(t)(r - \x\) 2 (6) 
the linear term vanishes, leaving 

72 (x) = (l + ^) (r - |x|) 2 + 0(m|x|- 2 ) (7) 
for some t between Ixl and r. 
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In dimensions d > 3, from Q we have 

Jd(x) = \x\ 2 + a d m\x\ 2 ~ d + 0(m\x\~ d ). 

Setting 4>(x) = \x\ 2 + adin\x\ 2 ~ d , the linear term in the Taylor expansion (|6 
of (j) about \x\ = r again vanishes, yielding 

ld {x) = (l + (d- l)(r/t) d ) (r - \x\) 2 + 0(m\x\- d ) 

for t between \x\ and r. Together with Q, this yields the following estimates 
in all dimensions d > 2. 

Lemma 2.1. Let 7^ 6e given by For a// 1 £ Z d ioe /iai>e 



7^)>(r--k|) 2 + 0(-^). (8) 



^n r, 



Lemma 2.2. Let 7^ 6e given by Then uniformly 

7d(x)=0(l), x 6 B r+ i - B r -i. 

The following lemma is useful for x near the origin, where the error term 
in (|8| blows up. 

Lemma 2.3. Let j d be given by Then for sufficiently large r, we have 

ld{x) > — , Vx G B r / 3 . 



Proof. Since 7d(x) — |x| 2 is super harmonic, it attains its minimum in B r / S 
at a point 2 on the boundary Thus for any x £ L? r /3 

7d(x) - \x\ 2 > -y d {z) - \z\ 2 , 



hence by Lemma 2.1 



r 2 



7d (x) > (2r/3) 2 -(r/3)^ + 0(l) > -. □ 



Lemmas 2.1 and 2.3 together imply the following. 

Lemma 2.4. Let 7^ 6e given by T/iere is a constant a depending only 
on d, such that j d > —a everywhere. 
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3 Divisible Sandpile 



Let no be a nonnegative function on Z rf with finite support. We start with 
mass Ho(y) at each site y. The operation of toppling a vertex x yields the 
mass distribution 

T x (i = Ho + a(x)A6 x 

where a(x) = max(fio(x) — 1,0) and A is the discrete Laplacian on Z d . Thus 
if Ho( x ) < 1 then T x ho = A*o an d no mass topples; if Ho( x ) > 1 then the 
mass in excess of 1 is distributed equally among the neighbors of x. 

Let x\, X2, ■ ■ ■ € Z d be a sequence with the property that for any x £ Z d 
there are infinitely many terms x k = x. Let 

Wfc(y) = T x k ■ ..T Xl Ho(y). 

be the amount of mass present at y after toppling the sites x±, . . . ,x k in 
succession. The total mass emitted from y during this process is 

My) : = Yl Vj-i(y) -Hj{v) = Yl ^ 

j<k:xj=y 3<k:x 3 =y 

where ctj(y) = m.ax(fij(y) — 1,0). 

Lemma 3.1. As k | oo the functions u k and Hk tend to limits u k ] u and 
Hk — ► H- Moreover, these limits satisfy 

H = Ho + < 1. 

Proof. Write M = ^ x Ho{ x ) f° r the total starting mass, and let B C Z d 
be a ball centered at the origin containing all points within L 1 -distance M 
of the support of ho- Note that if Hk{x) > and Ho( x ) = 0, then x must 
have received its mass from a neighbor, so Hk{y) > 1 for some y ~ x. Since 
*}2 x Hk{x) = M, it follows that ^ is supported on 1?. Let R be the radius 
of B, and consider the quadratic weight 

Qk=Yl Mx)\x\ 2 < MR 2 . 

x&L d 

Since Hk(x k )-Hk-i(x k ) = -a k (x k ) and for y ~ x k we have = 
^a k (x k ), we obtain 

Qfc - Qfc-i = «fc(Xfc) f — X 1 2/ 1 2 - kfcl 2 J = "fcOfc)- 
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Summing over k we obtain from ([9| 

Qk = Qo + u k( x )- 

Fixing x, the sequence Uk(x) is thus increasing and bounded above, hence 
convergent. 

Given neighboring vertices x ~ y, since y emits an equal amount of mass 
to each of its 2d neighbors, it emits mass Uk{y) /2d to x up to time k. Thus 
x receives a total mass of ^ ^Z y ^ x u k{y) from its neighbors up to time k. 
Comparing the amount of mass present at x before and after toppling, we 
obtain 

/ifc(x) = u (x) + Au k (x). 

Since j u we infer that fij. — > fj, := fj,o + Am. Note that if x^ = x, then 
Hk( x ) — 1- Since for each x£Z' ! this holds for infinitely many values of k, 
the limit satisfies \i < 1. □ 

A function s on Z d is superharmonic if As < 0. Given a function 7 on 
Z d the least superharmonic majorant of 7 is the function 

s(x) = inf{/(x) I / is superharmonic and / > 7}. 

The study of the least superharmonic majorant is a classical topic in analysis 
and PDE; see, for example, jS]. Note that if / is superharmonic and / > 7 
then 

y~x y~x 

Taking the infimum on the left side we obtain that s is superharmonic. 
Lemma 3.2. The limit u in Lemma \3.1\ is given by u = s + 7, where 

7(2) = \x\ 2 + y)na{y) 
and s is the least superharmonic majorant of —7. 



Proof. By Lemma 3.1 we have 

An = fj, — no < 1 — hq. 

Since A7 = 1 — no, the difference u— 7 is superharmonic. As u is nonnegative, 
it follows that n — 7 > s. For the reverse inequality, note that s + 7 — n is 
superharmonic on the domain D = {x \ fi(x) = 1} of fully occupied sites 
and is nonnegative outside D, hence nonnegative inside D as well. □ 
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As a corollary of Lemmas 3.1 and 3.2, we obtain the abelian property of 
the divisible sandpile, Proposition 1.2, which was stated in the introduction. 

We now turn to the case of a point source mass m started at the origin: 
/Uo = m5 . More general starting distributions are treated in [16], where 
we identify the scaling limit of the divisible sandpile model and show that 
it coincides with that of internal DLA and of the rotor-router model. In 
the case of a point source of mass m, the natural question is to identify 
the shape of the resulting domain D m of fully occupied sites, i.e. sites x for 
which n(x) = 1. According to Theorem 3.3, D m is extremely close to a ball 
of volume m; in fact, the error in the radius is a constant independent of m. 
As before, for r > we write 

B r = {xeZ d : |x| < r} 

for the lattice ball of radius r centered at the origin. 

Theorem 3.3. For m > let D m C Z d be the domain of fully occupied sites 
for the divisible sandpile formed from a pile of size m at the origin. There 
exist constants c, c' depending only on d, such that 

B r _ c C D m C B r+C i, 

where r = (m/uJd) 1 ^ and oJd is the volume of the unit ball in ~R d . 



The idea of the proof is to use Lemma 3.2 along with the basic estimates 



on 7, Lemmas |2.1| and 2.2 to obtain estimates on the odometer function 

u(x) = total mass emitted from x. 
We will need the following simple observation. 

Lemma 3.4. For every point x G D m — {o} there is a path x = xq ~ x\ ~ 
. . . ~ Xk = o in D m with n(xj + i) > u(xi) + 1. 

Proof. If Xi G D m — {o}, let Xi + \ be a neighbor of x; L maximizing 
Then G D m and 



u{x i+1 ) > — u (y) 



2d 



= u(xi) + Au(xi) 
= u(xi) + 1, 

where in the last step we have used the fact that x% G D m . □ 
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Proof of Theorem 3.3 We first treat the inner estimate. Let jd be given by 
By Lemma 3.2 the function u — 7^ is superharmonic, so its minimum in 



the ball B r is attained on the boundary. Since u > 0, we have by Lemma 2.2 
u(x) — 7rf(x) > — C, x G dB r 



for a constant C depending only on d. Hence by Lemma 2.1 
u(x) > (r - \x\) 2 - C'r d /\x\ d , x £ B r . 



(10) 



for a constant C depending only on d. It follows that there is a constant c, 
depending only on d, such that u{x) > whenever r/3 < \x\ < r — c. Thus 
-B r _ c — B r / 3 C D m . For 2; 6 -B r /3, by Lemma 2.3 we have u(x) > r 2 /4 — C > 
0, hence B r /^ C -D m . 



For the outer estimate, note that u— 7^ is harmonic on By Lemma 2.4 



we have 7^ > —a everywhere, where a depends only on d. Since u vanishes 
on dD m it follows that u — 7^ < a on D m . Now for any x £ D m with 
r — 1 < |x| < r, we have by Lemma |2.2| 

u(x) < 7d(x) + a < c' 



for a constant d depending only on d. Lemma 3.4 now implies that D m C 
Br+c'+l- n 



4 Classical Sandpile 

We consider a generalization of the classical abelian sandpile, proposed by 
Fey and Redig |5j. Each site in 7L d begins with a "hole" of depth H. Thus, 
each site absorbs the first H grains it receives, and thereafter functions 
normally, toppling once for each additional 2d grains it receives. If H is 
negative, we can interpret this as saying that every site starts with h = —H 
grains of sand already present. Aggregation is only well-defined in the regime 
h < 2d — 2, since for h = 2d — 1 the addition of a single grain already causes 
every site in Z d to topple infinitely often. 

Let S nt H be the set of sites that are visited if n particles start at the 
origin in Z d . Fey and Redig Theorem 4.7] prove that 

lim limsup — # {S n ,H A B H -i /dr ) = 0, 

where n = u>d rd , an d A denotes symmetric difference. The following theo- 
rem strengthens this result. 
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Theorem 4.1. Fix an integer H > 2 — 2d. Let S n = S n ^H be the set of sites 
that are visited by the classical abelian sandpile model in Z rf , starting from 
n particles at the origin, if every lattice site begins with a hole of depth H. 
Write n = uJd rd - Then 

where 

a = {2d-i + H)- 1 ' d 

and C2 is a constant depending only on d. Moreover if H > 1 — d, then for 
any e > we have 

S n ,H C #cir+4 

where 

c' 1 = (d-e + H)-V d 
and c' 2 is independent of n but may depend on d, H and e. 

Note that the ratio ci/c^ f 1 as H f do. Thus, the classical abelian 
sandpile run from an initial state in which each lattice site starts with a 
deep hole yields a shape very close to a ball. Intuitively, one can think of the 
classical sandpile with deep holes as approximating the divisible sandpile, 
whose limiting shape is a ball by Theorem 3.3 Following this intuition, 
we can adapt the proof of Theorem 3.3 to prove Theorem 4.1 just one 
additional averaging trick is needed, which we explain below. 

Consider the odometer function for the abelian sandpile 

u(x) = total number of grains emitted from x. 

Let T n = {x\u(x) > 0} be the set of sites which topple at least once. Then 

T n C S n C T n UdT n . 

In the final state, each site which has toppled retains between and 2d — 1 
grains, in addition to the H that it absorbed. Hence 

H < Au(x) + nb ox < 2d - 1 + H, x G T n . (11) 

We can improve the lower bound by averaging over a small box. For x € Z d 
let 

Qk(x) = {y£Z d : ||as — y||oo < k} 
be the box of side length 2k + 1 centered at x, and let 

u (k \x) = (2k + iy d <y)- 

y&Qk{x) 
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Write 

T h k) = i x I Qk{x) C T n }. 



Le Borgne and Rossin (12] observe that if T is a set of sites all of which 
topple, the number of grains remaining in T is at least the number of edges 
internal to T: indeed, for each internal edge, the endpoint that topples last 
sends the other a grain which never moves again. Since the box Qk(x) has 
2dk(2k + l)^" 1 internal ed ges, we have 



Au«(s) > ^-^d + H- ?7Tr ^l Qfc(o) (x), x e TW. (12) 



2k j rr n 

d + H - — — 

2k + 1 (2k + 1 

The following lemma is analogous to Lemma |3.4 



Lemma 4.2. For every point x E T n adjacent to dT n there is a path x = 
xq ~ x\ ~ . . . ~ x m = o in T n with u(xi+i) > u(xi) + 1. 

Proof. By (11) we have 

Since u(xj_i) < u(xi), some term u(y) in the sum above must exceed u(xi). 
Let x i+ i = y. □ 



Proof of Theorem 4-1 Let 

Ux) = {2d-l + H)\x\ 2 +ng{x), 

and let 

Cd{x) = £ d (x) - £d(lcirjei). 
Taking m = nf {2d — 1 + H) in Lemma |2.2| we have 

u(x)-Z d (x)>-U{x)>-C(2d-l + H), xe8B cir (13) 



for a constant C depending only on d. By (11), u — £,d is superharmonic, so 
u — £,d > — C(2d — 1 + H) in all of B Cir . Hence by Lemma 2.1 we have for 
x G B Cir 

u{x) >{2d-l + H) (W - |a;|) 2 - C \c x r) d /\xf) , (14) 

where C depends only on d. It follows that u is positive on B cir _ C2 — B Cir / s 
for a suitable constant C2 depending only on d. For x E B Cir / s , by Lemma 
we have u(sc) > (2d - 1 + H)(c\r 2 /A - C) > 0. Thus B Cir - C2 C T n C S n 



2.3 
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For the outer estimate, let 

4>d(x) = (d — e + H)\x\ 2 + ng(x). 
Choose k large enough so that ji^tid > d — e, and define 

Mx) = (2k + l)- d £ ^(y). 

2/eQ fe (a;) 

Finally, let 

^d(x) = Tpd(x) - Tpd([c'ir\ei). 
By (12), — ^ is subharmonic on T„ . Taking m = n/(d — e + -ff) 



in Lemma |2.4| there is a constant a depending only on d, such that ^ > 
-a(d + H) everywhere. Since < (2d + H) ( - d+l ^ k on <9Ti fc) it follows that 
- < o(d + H) + (2d + Hf d+1 ) k on r^ fc) . Now for any x G 5„ with 



c^r — 1 < |x| < d^r we have by Lemma 



2.2 



(fe) 



(z) < ^ d (x) + a(d + H) + (2d + F)( d+1 ) fc < c 2 



for a constant c 2 d epending only on d, H and e. Then u(z) < c' 2 := (2/c + 
l) d C2- Lemma 4.2 now implies that T n C i? c ' ir+c ^, and hence 



We remark that the crude bound of (2d + H)( d+1 ^ k used in the proof 
of the outer estimate can be improved to a bound of order k 2 H, and the 
final factor of (2k + l) d can be replaced by a constant factor independent 
of k and H, using the fact that a nonnegative function on Z rf with bounded 
Laplacian cannot grow faster than quadratically; see [16] . 

5 Rotor-Router Model 

Given a function / on Z d , for a directed edge (x, y) write 

Vf(x,y) = f(y)-f(x). 

Given a function s on directed edges in write 

divs(x) = — ^2s(x,y). 

The discrete Laplacian of / is then given by 

Af(x) = dwvf = 1 ^ d Yl f(y) - f^- 
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5.1 Inner Estimate 

Fixing n > 1, consider the odometer function for rotor-router aggregation 

u(x) = total number exits from x by the first n particles. 

We learned the idea of using the odometer function to study the rotor- router 
shape from Matt Cook [2 . 

Lemma 5.1. For a directed edge (x,y) in Z d , denote by n{x,y) the net 
number of crossings from x to y performed by the first n particles in rotor- 
router aggregation. Then 

Vu(x, y) = — 2dK,(x, y) + R(x, y) (15) 

for some edge function R which satisfies 

\R(x,y)\ <4d-2 

for all edges (x, y). 

Remark. In the more general setting of rotor stacks of bounded discrepancy, 
the 4d — 2 will be replaced by a different constant here. 

Proof. Writing N(x, y) for the number of particles routed from x to y, we 
have 

<x)-2d + l u(x) + 2d-l 

2d - Jv ^,yj- 2d 

hence 

\Vu(x, y) + 2dK(x, y)\ = \u{y) - u{x) + 2dN(x, y) - 2dN(y, x)\ 
< 4d - 2. □ 



In what follows, Co,C\, . . . denote constants depending only on d. 
Lemma 5.2. Let Q C Z d - {o} with 2<#n<oo. Then 

yen 
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Proof. For each positive integer k, let 

S k = {yeZ d :k< \y\ < k + 1}. 

Then 

y&S k 

for a constant C depending only on d. Since 0, can intersect at most 
Diam(S7) + 1 < 2Diam(S7) distinct sets <S&, taking Co = 2C the proof is 
complete. □ 

Lemma 5.3. Let G = Gb t be the Green's function for simple random walk 
in Ij d stopped on exiting B r . For any p > 1 and x £ B r , 

E I^lG^-G^I^Cip. (16) 

|z-y|<p 

Proof. Let (Xj)j>o denote simple random walk in Z rf , and let T be the first 
exit time from B r . For fixed y, the function 

A(x) = g(x -y)- E x g(X T - y) (17) 

has Laplacian AA(x) = —5 xy in B r and vanishes on dB r , hence A(x) = 
G(x,y). 

Let x,y £ B r and z ~ y. From Q we have 

3(z - y) - g(x - z) \< r^— f , y,z^x. 

\x — y\ a 1 



Using the triangle inequality together with (17), we obtain 
\G(x,y) - G(x, z)\ < \g(x -y)-g{x-z)\+ E x \g{X T - y) - g(X T - z)\ 

<^^+ E wAi- 

a; — y z — ' \w — y\ 

1 Wl wedB r 1 y| 

where H x (w) = P x (Xr = if). 

Write D = {y £ B r : \x — y\ <}. Then 

E£|G(x )2 /)-G(x,z)| <C 3 p + C 2 £ ^HE^-^I 1 ^ ( 18 ) 

yeD z ~V w£dB r y£D 



18 



Hi 



Hi 









\ 



Figure 4: Diagram for the Proof of Lemma 5.4 



Taking $7 = w — D in Lemma 5.2 the inner sum on the right is at most 



CoDiam(D) < 2Cop, so the right side of (18) is bounded above by C\p for 
a suitable C\. 

Finally, the terms in which y or z coincides with x make a negligible 
contribution to the sum in (16), since for y ~ x £ Z rf 



|G(a;,x)-G(x,y)| < | 5 (o) - s(x-y)| +E x \g(X T -x) - g(X 7 



< C 4 . □ 



Lemma 5.4. Lei i^i, i?2 be linear half-spaces in Z d , not necessarily parallel 
to the coordinate axes. Let Ti be the first hitting time of Hi. If x ^ H\ UH2, 
then 

™ /m m \ 5/11+I / 1 

> T 2 ) < 



2 /i 2 

where hi is the distance from x to Hi. 



1 + 



2hn 



Proof. If one of H\, H2 contains the other, the result is vacuous. Otherwise, 
let Hi be the half-space shifted parallel to Hf by distance 2ti2 in the direction 
of x, and let Tj be the first hitting time of HiUHi. Let (Xt)t>o denote simple 
random walk in Z rf , and write Mt for the (signed) distance from Xt to the 
hyperplane defining the boundary of Hi, with Mo = hi. Then Mt is a 
martingale with bounded increments. Since E^Ti < 00, we obtain from 
optional stopping 



hi = E x M fi > 2h 2 F x {X fi £ H{) - F»(X fi £ Hi), 
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hence 



(X? G #i) < 



hi + 1 

2h 2 



(19) 



Likewise, dM 2 — t is a martingale with bounded increments, giving 



1 1 



< d{2h 2 + 1) 2 {x fi e H x ) 



<d(h 1 + l)(2h 2 + l) 1 + 



2/i 2 



(20) 



Let T = min(Ti,T2). Denoting by Z?i the distance from X% to the 
hyperplane defining the boundary of H 2 , the quantity 

N t = ^{D 2 + {2h 2 -D t ) 2 )-t 

is a martingale. Writing p = P X (T = T 2 ) we have 

dh\ = EiV = EiVr > P^(2h 2 ) 2 + (1 - p)d/i| - E^T 
> (1 + p)dh 2 — E X T 



hence by ( 20 ) 



P< 



E^T 
dh\ 



< 2 



hi + 1 



1 + 



2/i 5 



Finally by (19) 



P(T 1 >T 2 )< P + P(x fi efli)< 



5 hi + 1 

2 /l 2 



1 + 



2/l 5 



□ 



Lemma 5.5. Let x £ B r and let p = r + 1 — \x\. Let 

S* k = {ye B r : 2 k p <\x-y\< 2 k+l p}. (21) 
Let Tfc be the first hitting time of St, and T the first exit time from B r . Then 

W x (T k <T)< C 2 2- k . 
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Figure 5: Diagram for the proof of Lemma 



5.5 



Proof. Let H be the outer half-space tangent to B T at the point z £ dB r 
closest to x. Let Q be the cube of side length 2 k ' pj\fd centered at x. Then 
Q is disjoint from Sj*, hence 

W x (T k <T)< F x (T dQ <T)< F x (T dQ < T H ) 

where Tqq and Tjj are the first hitting times of dQ and H. Let Hi, . . . , i?2d 
be the outer half-spaces defining the faces of Q, so that Q = Hf n . . . n 
By Lemma 5.4 we have 

2d 

W x {T dQ < T H ) < Y, F ^ T H t < T H ) 
i=l 

„ 5 ^ dist(x, H) + 1 / 1 



2^ dist(x,Hi) V 2dist(x,ff 1 



Since dist(x,-ff) = \x — z\ < p and dist(a?, Hi) = 2 k 1 p/\/ r d, and p > 1, 
taking C 2 = 20 <i 3 / 2 (l + \/d) 2 completes the proof. □ 

Lemma 5.6. Let G = Gb t be the Green's function for random walk stopped 
on exiting B r . Let x £ B r and let p = r + 1 — |x| . Then 

Y, ^|G(x,y)-G(x,z)| <C 3 plog-. 

yeB r z~y P 
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Proof. Let S k be given by (21 ), and let 

W = {w e d(S* k UdS* k ) : \w-x\<2 k p} 

be the portion of the boundary of the enlarged spherical shell S k U dS k lying 
closer to x. Let T\\r be the first hitting time of W, and T the first exit time 
from B r . For w S W let 

H x {w)=F x {X TwAT = w). 

For any y E <S£ and 2 ~ y, simple random walk started at x must hit W 
before hitting either y or z, hence 

\G(x,y)-G(x,z)\ < H x (w)\G(w,y)-G(w,z)\. 

For any y £ S k and any w £ W we have 

\y — w\ < \y — x\ + \ w — x\ < 3 ■ 2 k p. 

Lemma |5.3| yields 

X)|G(x,y) -<?(*,*) | < 3d2 fc p ^ i^H- 



By Lemma 5.5 we have Ylwew Hx{ w ) < C-22 , so the above sum is at most 
3C1C2/). Since the union of shells Sq , iSj*, . . . , <5*i og ( r / p y\ covers all of B r ex- 
cept for those points y within distance p of x, and Yl\ y - X \<p Yl z ~ y \G(x, y) — 
G(x, z)\ < C\p by Lemma 5.3, the result follows. □ 

Proof of Theorem Inner Estimate. Let n and R be defined as in Lemma 



5.1 Since the net number of particles to enter a site x 7^ o is at most one, we 



have 2d div k(x) > —1. Likewise 2d div k(o) = n — 1. Taking the divergence 
in (15 ), we obtain 

Au(x) < l + divi?(x), x 7^0; (22) 
Au(o) = 1 - n + div i?(o) . (23) 

Let T be the first exit time from B r , and define 

f{x) = E x u{X T ) — E X T + n E x #{j < T\Xj = 0}. 

Then A/(x) = 1 for x € B r — {0} and A/(o) = 1 — n. Moreover / > on 
dB r . It follows from Lemma |2.2| with m = n that / > 7 — C4 on 5 r for a 
suitable constant C4. 
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We have 



u{x) - E x u(X T ) = ^2E x (u(X kAT ) - u(X( k+1 



k>0 



Each summand on the right side is zero on the event {T < k}, hence 
E x (u(X kAT ) -n(X (fc+1)AT ) | F k AT) = -Au(X k )l {T>k} . 



Taking expectations and using (|22|) and (23), we obtain 



u(x) - E x u{X T ) > J2 E x[l{T>k}{nl {Xk =o} ~ l-divR{X k )) 



k>0 

n E x #{k <T\X k = o}- E X T -^E x [l {T>fc} div R(X k ) 

k>0 



hence 



k>0 



z~X k 



(24) 



Since random walk exits B r with probability at least ^ every time it reaches 
a site adjacent to the boundary dB r , the expected time spent adjacent to 



the boundary before time T is at most 2d. Since \R\ < 4d, the terms in (24) 
with z £ dB r contribute at most 16c? 3 to the sum. Thus 



k>0 



1 {T>k}n{X h =y}R(y,z) 



y,z£B r 
y~z 



8d 2 . 



For y £ B r we have {X k = y} D {T > k] = {X kA x = y}, hence 

u(x)-f(x)>-^Y, E V x (X kAT = y)R(y,z)-8d 2 . (25) 



k>0 y,z£B r 
y^z 



Write p k {y) = E x (X kA j> = y). Note that since Vii and n are antisymmet- 
ric, R is antisymmetric. Thus 



y,zGB r 
y~z 



y,z(LB r 



y,z(^B r 
y~z 
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Summing over k and using the fact that \R\ < Ad, we conclude from (25) 
that 

u(x)>f(x)- \G(x,y)-G(x,z)\-8d 2 , 

y,z&B r 

where G = Gb t is the Green's function for simple random walk stopped on 
exiting B r . By Lemma |5,6| we obtain 

u(x) > f{x) - C 3 (r + 1 - \x\) log , /" , , - 8a! 2 . 

r + 1 — \x\ 



Using the fact that / > 7 — C4, we obtain from Lemma 2.1 



u(x) > (r - \x\) 2 - C 3 (r + 1 - |x|) log r + O ( ^~ 

r + 1 — \x\ \ \x\ a 

The right side is positive provided r/3 < \x\ < r — C^logr. For x 6 -B r /3, by 
we have u(x) > r 2 /A — C^r log | > 0, hence B r _c 5 \ogr C ^4 n - D 



Lemma 



2.3 



5.2 Outer Estimate 

The following result is due to Holroyd and Propp (unpublished); we include 



a proof for the sake of completeness. Notice that the bound in (26) does not 
depend on the number of particles. 

Proposition 5.7. Let T be a finite connected graph, and letYcZ be sub- 
sets of the vertex set ofT. Let s be a nonnegative integer-valued function on 
the vertices ofT. Let H w (s,Y) be the expected number of particles stopping 
in Y if s(x) particles start at each vertex x and perform independent sim- 
ple random walks stopped on first hitting Z. Let H r (s,Y) be the number of 
particles stopping in Y if s(x) particles start at each vertex x and perform 
rotor-router walks stopped on first hitting Z. Let H{x) = H W (1 X ,Y). Then 

\H r (s,Y) - H w (s,Y)\ <J2J2\ H W- H ( V> >\ ( 26 ) 

independent of s and the initial positions of the rotors. 

Proof. For each vertex u ^ Z, arbitrarily choose a neighbor rj(u). Or- 
der the neighbors r](u) = v 1, t>2, . . . , of w so that the rotor at u points 
to Vi + i immediately after pointing to V{ (indices mod d). We assign weight 
w(u,r](u)) = to a rotor pointing from u to r)(u), and weight w(u,Vi) = 
H(u) — H(vi) + w(u,Vi-i) to a rotor pointing from u to V{. These assign- 
ments are consistent since H is a harmonic function: Y2i(H( u ) ~^i v i)) = 0- 
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Figure 6: Diagram for the proof of Lemma 



5.8 



We also assign weight H(u) to a particle located at u. The sum of rotor 
and particle weights in any configuration is invariant under the operation of 
routing a particle and rotating the corresponding rotor. Initially, the sum of 
all particle weights is H w (s, Y). After all particles have stopped, the sum of 
the particle weights is H r (s, Y). Their difference is thus at most the change 



in rotor weights, which is bounded above by the sum in (26). □ 
For peZ let 

S p = {x £ Z d : p < \x\ < p + 1}. (27) 

Then 

B p = {xeZ d : \x\ <p} = 5 U...U5p_i. 

Note that for simple random walk started in B p , the first exit time of B p 
and first hitting time of S p coincide. Our next result is a modification of 
Lemma 5(b) of |10j . 



Lemma 5.8. Fix p > 1 and y G S p . For x € B p let H(x) = F x (Xt = y), 
where T is the first hitting time of S p . Then 

H{X) ~ \x-y\^ (28) 
for a constant J depending only on d. 
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Proof. We induct on the distance \x — y\, assuming the result holds for all 
x' with \x' — y\ < 7}\x — y\; the base case can be made trivial by choosing J 
sufficiently large. By Lemma 5(b) of (IQj . we can choose J large enough so 
that the result holds provided |y| — \x\ > 2~ d ~ 3 \x — y\. Otherwise, let H\ 
be the outer half-space tangent to S p at the point of S p closest to x, and let 
H2 be the inner half-space tangent to the ball S of radius — y\ about y, 



at the point of S closest to x. By Lemma 5.4 applied to these half-spaces, 
the probability that random walk started at x reaches S before hitting S p 
is at most 2 l ~ d . Writing T for the first hitting time of S U S p , we have 



H(x) < W x {X f = x')H(x') < 2 l - d J 



\x - y\ 



l-d 



x'es 

where we have used the inductive hypothesis to bound H{x'). □ 

The lazy random walk in Z d stays in place with probability ^, and moves 
to each of the 2d neighbors with probability • We will need the following 
standard result, which can be derived e.g. from the estimates in [IT], section 
11.12; we include a proof for the sake of completeness. 

Lemma 5.9. Given u ~ v G Z , lazy random walks started at u and v can 
be coupled with probability 1 — C/R before either reaches distance R from u, 
where C depends only on d. 

Proof. Let i be the coordinate such that m 7^ V{. To define a step of the 
coupling, choose one of the d coordinates uniformly at random. If the chosen 
coordinate is different from i, let the two walks take the same lazy step so 
that they still agree in this coordinate. If the chosen coordinate is i, let 
one walk take a step while the other stays in place. With probability \ the 
walks will then be coupled. Otherwise, they are located at points u', v' with 
\u' — v'\ = 2. Moreover, P(|n — v!\ > < % for a constant C depending 
only on d. From now on, whenever coordinate i is chosen, let the two walks 
take lazy steps in opposite directions. 
Let 

_ u[ + v \ 
<~ 2 



be the hyperplane bisecting the segment [u' , v']. Since the steps of one walk 
are reflections in H\ of the steps of the other, the walks couple when they 
hit H\. Let Q be the cube of side length R/^/d + 2 centered at u, and let H2 



be a hyperplane defining one of the faces of Q. By Lemma 5.4 with h\ = 1 
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and /12 = R/4y/d, the probability that one of the walks exits Q before the 



walks couple is at most 2d 



5 hi+1 
2 h 2 



< 40d 3 / 2 (l + 2Vd) 2 /R. 



□ 



Lemma 5.10. With H defined as in Lemma 5.8 we have 

H(v)\ < J'logp 



U 



for a constant J' depending only on d. 



Proof. Given u 6 B p and v ~ u, by Lemma 5.9 lazy random walks started 
at u and v can be coupled with probability 1 — 2C/\u — y\ before either 
reaches distance \u — y\/ 2 from u. If the walks reach this distance without 
coupling, by Lemma 5.8 each has still has probability at most J/\u — y\ d ~ l 
of exiting B p at y. By the strong Markov property it follows that 

2C T 

\H(u) - H(v)\ < v ^ J ld . 

\u — y\ a 

Summing in spherical shells about y, we obtain 

ueBp v~u t=l 



We remark that Lemma 15.101 could also be inferred from Lemma 15.81 
using [9, Thm. 1.7.1] in a ball of radius \u — y\/2 about u. 

To prove the outer estimate of Theorem |1.1| we will make use of the 
abelian property of rotor-router aggregation. Fix a finite set F C con- 
taining the origin. Starting with n particles at the origin, at each time step, 
choose a site x G F with more than one particle, rotate the rotor at x, and 
move one particle from x to the neighbor the rotor points to. After a finite 
number of such choices, each site in F will have at most one particle, and all 
particles that exited F will be on the boundary dF. The abelian property 
says that the final configuration of particles and the final configuration of 
rotors do not depend on the choices. For a proof, see (3J Prop. 4.1]. 

In our application, we will fix p > r and stop each particle in rotor- router 
aggregation either when it reaches an unoccupied site or when it reaches the 
spherical shell S p . Let N p be the number of particles that reach S p during 
this process. Note that at some sites in S p , more than one particle may have 
stopped. If we let each of these extra particles in turn continue performing 
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rotor-router walk, stopping either when it reaches an unoccupied site or 
when it hits the larger shell S p+ h, then by the abelian property, the number 
of particles that reach S p +h will be N p+ h. We will show that when h is order 
r 1-1 / d , a constant fraction of the particles that reach S p find unoccupied sites 
before reaching S p+ h- 

Proof of Theorem \l.l\ Outer Estimate. Fix integers p > r and h > 1. In 
the setting of Proposition 5.7, let T be the lattice ball B p+ i %+ i, and let 
Z = S p+ h- Fix y £ S p+ h and let Y = {y}. For x £ S p , let s(x) be the 
number of particles stopped at x if each particle in rotor-router aggregation 
is stopped either when it reaches an unoccupied site or when it reaches S p . 
Write 

H(x) = F X (X T = y) 



where T is the first hitting time of S p+ h- By Lemma 5.8 we have 



H w (s,y) = s(x)H(x) < 
xeSp 



JNp 
h d ~ l 



(29) 



where 



is the number of particles that ever visit the shell S p . 
By Lemma 5.10 the sum in (26) is at most J'lo. 

ton 



5.7 and (29) we have 



By Lemma 5.10 the sum in (26) is at most J'log/i, hence from Proposi- 

(30) 



H r (s,y)<(^ + J'logh. 



h d ~ 

Let p(0) = r, and define p{i) inductively by 

p{i + 1) = min [p{i) + N%g d - 1) , min{p > p(i)\N p < N p(i) /2}} . (31) 
Fixing h < p{i + 1) — p(i), we have 

2d-2 

h d - l \ogh< N^ 1 logN p{l) < N p{l y, 



so (30) with p = p{i) simplifies to 



H r (s,y) < 



CN, 



P(i) 



(32) 



where C = J + J'. 



2S 



Since all particles that visit S p ^ + h during rotor- router aggregation must 
pass through S p u-\, we have by the abelian property 



N p(l)+h < Yl H r(^y)- 



(33) 



Let Mfc = #(A n f}Sk). There are at most M p ^ +h nonzero terms in the sum 
on the right side of ( 33 ) , and each term is bounded above by ( 32 ) , hence 

h d ~ l h d ~ l 
M p(t)+h > JV^^- > — 

where the second inequality follows from N p ^ +h > N p u\/2. Summing over 
h, we obtain 



p= P (i)+i 



i 

p ~ 2dC 



+ 1) _ p (i) - if. 



(34) 



The left side is at most N p ^, hence 



p(i + 1) - p(<) < {2dCN p[i) f/ d < N 2 p $ d -V 

provided N p{i) > C := (2dC) 2d ~ 1 . Thus the minimum in d3l| is not attained 
by its first argument. It follows that N p u + i\ < N p u\/2, hence N p ^ a i ogr ^ < C 
for a sufficiently large constant a. 

By the inner estimate, since the ball B r _ c \ ogr is entirely occupied, we 
have 

Yj M p ^ u dr d - Mr ~ clogr) d 

p>r 

< cdujdr^ 1 log r. 



Write Xi = p{i + 1) — p{i) — 1; by (34) we have 

a log r 

V] xf < cduJdr^ 1 logr, 



i=0 



By Jensen's inequality, subject to this constraint, ^ Xi is maximized when 
all equal, in which case Xi < C"r 1 - x / d and 



Since N 



p(alogr) 



p(alogr) = r + X{ < r + C"r l 1 ^ d logr. (35) 
< C we have N p r a i ogr \ + o> = 0; that is, no particles reach 



the shell 5 p ( a log r)+c • Taking d = C + C", we obtain from (35) 

A n C -B r (l+ c ' r -i/di ogr .). □ 



29 



Figure 7: Image of the rotor-router aggregate of one million particles under 
the map z ^ 1/z 2 . The colors represent the rotor directions. The white 
disc in the center is the image of the complement of the occupied region. 



6 Concluding Remarks 

A number of intriguing questions remain unanswered. Although we have 
shown that the asymptotic shape of the rotor-router model is a ball, the 
near perfect circularity found in Figure [T] remains a mystery. In particular, 



we do not know whether an analogue of Theorem 1.3 holds for the rotor- 
router model, with constant error in the radius as the number of particles 
grows. 

Equally mysterious are the patterns in the rotor directions evident in 
Figure [T] The rotor directions can be viewed as values of the odometer 
function mod 2d, but our control of the odometer is not fine enough to 
provide useful information about the patterns. If the rescaled occupied 
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region yTrfn A n is viewed as a subset of the complex plane, it appears that 
the monochromatic regions visible in Figure [T] in which all rotors point 
in the same direction, occur near points of the form (1 + 2z) -1 / 2 , where 
z = a + bi is a Gaussian integer (i.e. a, b £ Z). We do not even have 
a heuristic explanation for this phenomenon. Figure 7 shows the image 
of ^1,000,000 under the map z i— > 1/z 2 ; the monochromatic patches in the 
transformed region occur at lattice points. 

Laszlo Lovasz (personal communication) has asked whether the occu- 
pied region A n is simply connected, i.e. whether its complement is con- 
nected. While Theorem |1.1| shows that A n cannot have any holes far from 
the boundary, we cannot answer his question at present. 

A final question is whether our methods could be adapted to internal 
DLA to show that if n = uj^r , then with high probability -B r _ c i gr C I n , 
where I n is the internal DLA cluster of n particles. The current best bound is 
due to Lawler who proves that with high probability -6 r _ r i/3(i g r .)2 
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